Poisson binomial distribution (poisson_binom)#
The Poisson binomial distribution is the distribution of the sum of independent Bernoulli random variables with different success probabilities.
If \(X_i \sim \mathrm{Bernoulli}(p_i)\) independently for \(i=1,\dots,n\), then
Learning goals#
Recognize Poisson binomial data and typical modeling patterns.
Write the PMF and CDF, and understand the generating-function view.
Compute mean/variance/skewness/kurtosis via additivity of cumulants.
Implement PMF/CDF and sampling (NumPy-only) and visualize behavior.
Use
scipy.stats.poisson_binomfor probability calculations and hypothesis tests.
Prerequisites#
Expectation and variance (linearity and independence)
Comfort with logs, products, and basic numerical computing
Notebook roadmap#
Title & Classification
Intuition & Motivation
Formal Definition
Moments & Properties
Parameter Interpretation
Derivations (Expectation, Variance, Likelihood)
Sampling & Simulation (NumPy-only)
Visualization (PMF, CDF, Monte Carlo)
SciPy Integration
Statistical Use Cases
Pitfalls
Summary
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
from scipy import stats
from scipy.optimize import minimize
from scipy.special import expit, xlogy
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=4, suppress=True)
SEED = 7
rng = np.random.default_rng(SEED)
import sys
import scipy
import plotly
print("Python:", sys.version.split()[0])
print("NumPy:", np.__version__)
print("SciPy:", scipy.__version__)
print("Plotly:", plotly.__version__)
print("Seed:", SEED)
Python: 3.12.9
NumPy: 1.26.2
SciPy: 1.15.0
Plotly: 6.5.2
Seed: 7
1) Title & Classification#
Name:
poisson_binom(Poisson binomial distribution)Type: Discrete
Support: \(k \in \{0,1,\dots,n\}\) where \(n = \lvert\mathbf p\rvert\) is the number of Bernoulli trials
Parameter space:
\(\mathbf p = (p_1,\dots,p_n)\) with \(0 \le p_i \le 1\) for all \(i\).
(SciPy) an integer shift
locis also available; it shifts the support to \(\{\mathrm{loc}, \dots, \mathrm{loc}+n\}\).
Notation:
\(X \sim \mathrm{PoissonBinomial}(p_1,\dots,p_n)\).
Sometimes written as \(X \sim \mathrm{PB}(\mathbf p)\).
2) Intuition & Motivation#
What this distribution models#
You observe a count of successes across \(n\) Bernoulli trials, but the trials are heterogeneous:
Trial \(i\) succeeds with probability \(p_i\).
Trials are independent, but the \(p_i\) are not necessarily equal.
So \(X\) is a heterogeneous version of a binomial count.
Typical real-world use cases#
Reliability engineering: number of component failures when components have different failure probabilities.
Marketing / product analytics: number of conversions when each user/session has a different conversion probability.
Elections / forecasting: number of wins across districts with different win probabilities.
Multiple testing: number of discoveries when each test has a different Type-I error probability (or power).
Relations to other distributions#
Binomial: if all \(p_i\) are equal to a common \(p\), then \(X \sim \mathrm{Binomial}(n,p)\).
Poisson approximation: if all \(p_i\) are small and \(\lambda = \sum_i p_i\) is moderate, then \(X \approx \mathrm{Poisson}(\lambda)\).
Normal approximation (CLT): for large \(n\) with non-degenerate variance, \(X\) is approximately normal with mean \(\sum_i p_i\) and variance \(\sum_i p_i(1-p_i)\).
Sum of independent Bernoulli: conceptually, Poisson binomial is the distribution of that sum.
3) Formal Definition#
Let \(X_1,\dots,X_n\) be independent with
Define
PMF#
For \(k \in \{0,1,\dots,n\}\),
A very useful equivalent view is via the probability generating function (PGF):
Then \(\mathbb P(X=k)\) is exactly the coefficient of \(z^k\) in \(G_X(z)\):
CDF#
The CDF is the partial sum of the PMF:
There is no single closed form for general \((p_1,\dots,p_n)\), but it is efficiently computable.
def validate_p_vector(p) -> np.ndarray:
p = np.asarray(p, dtype=float)
if p.ndim != 1:
raise ValueError("p must be a 1D array-like of probabilities")
if p.size == 0:
raise ValueError("p must have at least one element")
if not np.all(np.isfinite(p)):
raise ValueError("p must be finite")
if np.any((p < 0.0) | (p > 1.0)):
raise ValueError("each p_i must be in [0, 1]")
return p
def poisson_binom_support(p: np.ndarray) -> np.ndarray:
p = validate_p_vector(p)
return np.arange(p.size + 1)
def poisson_binom_pmf_array(p: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Return (k, pmf[k]) for k=0..n using a stable O(n^2) recursion.
Recurrence: after processing probabilities p_1..p_i,
P_i(k) = P_{i-1}(k) (1-p_i) + P_{i-1}(k-1) p_i.
"""
p = validate_p_vector(p)
n = p.size
pmf = np.zeros(n + 1, dtype=float)
pmf[0] = 1.0
for pi in p:
pmf_next = np.zeros_like(pmf)
pmf_next[0] = pmf[0] * (1.0 - pi)
pmf_next[1:] = pmf[1:] * (1.0 - pi) + pmf[:-1] * pi
pmf = pmf_next
pmf = np.clip(pmf, 0.0, 1.0)
pmf = pmf / pmf.sum() # guard against tiny rounding drift
k = np.arange(n + 1)
return k, pmf
def poisson_binom_logpmf(k, p: np.ndarray) -> np.ndarray:
p = validate_p_vector(p)
n = p.size
k_arr = np.asarray(k)
out = np.full(k_arr.shape, -np.inf, dtype=float)
k_int = k_arr.astype(int)
valid = (k_int == k_arr) & (k_int >= 0) & (k_int <= n)
if not np.any(valid):
return out
_, pmf = poisson_binom_pmf_array(p)
out[valid] = np.log(pmf[k_int[valid]])
return out
def poisson_binom_pmf(k, p: np.ndarray) -> np.ndarray:
return np.exp(poisson_binom_logpmf(k, p))
def poisson_binom_cdf(x, p: np.ndarray) -> np.ndarray:
p = validate_p_vector(p)
n = p.size
x_arr = np.asarray(x)
out = np.zeros_like(x_arr, dtype=float)
out[x_arr >= n] = 1.0
inside = (x_arr >= 0) & (x_arr < n)
if np.any(inside):
k = np.floor(x_arr[inside]).astype(int)
_, pmf = poisson_binom_pmf_array(p)
cdf = np.cumsum(pmf)
out[inside] = cdf[k]
return out
4) Moments & Properties#
Because \(X\) is a sum of independent Bernoulli variables, many properties follow from additivity.
Mean and variance#
Using linearity of expectation and independence:
Higher standardized moments#
A convenient way to get skewness/kurtosis is via cumulants. For a Bernoulli\((p)\) variable, the first four cumulants are:
For independent sums, cumulants add: \(\kappa_r(X) = \sum_i \kappa_r(X_i)\).
Then
MGF and characteristic function#
The moment generating function (MGF) and characteristic function factorize:
Entropy#
There is no simple closed form in general. Given the PMF, entropy (in nats) is
def poisson_binom_cumulants(p: np.ndarray) -> dict:
p = validate_p_vector(p)
k1 = float(np.sum(p))
k2 = float(np.sum(p * (1.0 - p)))
k3 = float(np.sum(p * (1.0 - p) * (1.0 - 2.0 * p)))
k4 = float(np.sum(p * (1.0 - p) * (1.0 - 6.0 * p + 6.0 * p**2)))
return {"k1": k1, "k2": k2, "k3": k3, "k4": k4}
def poisson_binom_moments(p: np.ndarray) -> dict:
ks = poisson_binom_cumulants(p)
mean = ks["k1"]
var = ks["k2"]
if var == 0.0:
skew = float("nan")
excess_kurt = float("nan")
else:
skew = ks["k3"] / (var ** 1.5)
excess_kurt = ks["k4"] / (var ** 2)
return {
"mean": mean,
"var": var,
"skew": skew,
"excess_kurt": excess_kurt,
}
def poisson_binom_log_mgf(t, p: np.ndarray) -> np.ndarray:
"""log M_X(t) = sum log((1-p_i) + p_i e^t)."""
p = validate_p_vector(p)
t = np.asarray(t, dtype=float)
# logaddexp(log(1-p), log(p)+t) is stable at p=0 or p=1
log_terms = np.logaddexp(np.log1p(-p), np.log(p) + t[..., None])
return np.sum(log_terms, axis=-1)
def poisson_binom_mgf(t, p: np.ndarray) -> np.ndarray:
return np.exp(poisson_binom_log_mgf(t, p))
def poisson_binom_cf(omega, p: np.ndarray) -> np.ndarray:
p = validate_p_vector(p)
omega = np.asarray(omega)
return np.prod(1.0 - p + p * np.exp(1j * omega[..., None]), axis=-1)
def poisson_binom_entropy(p: np.ndarray, *, base=np.e) -> float:
p = validate_p_vector(p)
_, pmf = poisson_binom_pmf_array(p)
H_nats = -float(np.sum(xlogy(pmf, pmf)))
if base == np.e:
return H_nats
return H_nats / float(np.log(base))
p = np.array([0.10, 0.25, 0.60, 0.80])
k, pmf = poisson_binom_pmf_array(p)
cdf = np.cumsum(pmf)
mom = poisson_binom_moments(p)
print("p:", p)
print("support:", k)
print("pmf sum:", pmf.sum())
print("mean/var/skew/excess_kurt (from cumulants):", mom)
print("entropy (nats):", poisson_binom_entropy(p))
# Monte Carlo sanity check
samples = (rng.random((200_000, p.size)) < p).sum(axis=1)
print("MC mean:", samples.mean(), " | formula:", mom["mean"])
print("MC var :", samples.var(ddof=0), " | formula:", mom["var"])
p: [0.1 0.25 0.6 0.8 ]
support: [0 1 2 3 4]
pmf sum: 1.0
mean/var/skew/excess_kurt (from cumulants): {'mean': 1.75, 'var': 0.6775, 'skew': 0.03900275742662621, 'excess_kurt': -0.17698560749445102}
entropy (nats): 1.2220616806986953
MC mean: 1.752005 | formula: 1.75
MC var : 0.6768534799749998 | formula: 0.6775
5) Parameter Interpretation#
The parameter vector \(\mathbf p=(p_1,\dots,p_n)\) is the list of per-trial success probabilities.
Each \(p_i\) controls how likely trial \(i\) is to contribute a 1 to the sum.
The mean depends only on the sum \(\sum_i p_i\).
The variance depends on \(\sum_i p_i(1-p_i)\).
For a fixed mean, variance is maximized when the \(p_i\) are all equal (a consequence of concavity of \(p(1-p)\)).
Making probabilities more extreme (closer to 0 or 1) typically reduces variance and makes the distribution more concentrated.
Because the parameters are a vector, many different shapes are possible even with the same mean.
def show_pmf_comparison(ps: dict[str, np.ndarray], *, title: str):
fig = go.Figure()
for name, pvec in ps.items():
k, pmf = poisson_binom_pmf_array(pvec)
fig.add_trace(go.Scatter(x=k, y=pmf, mode="lines+markers", name=name))
fig.update_layout(
title=title,
xaxis_title="k (number of successes)",
yaxis_title="P(X = k)",
legend_title="parameter set",
)
fig.show()
n = 30
mean_target = 0.30
p_equal = np.full(n, mean_target)
# Same mean, more heterogeneous probabilities
p_mixture = np.r_[np.full(n // 2, 0.10), np.full(n - n // 2, 0.50)] # mean = 0.30
p_extreme = np.r_[np.zeros(n // 2), np.full(n - n // 2, 0.60)] # mean = 0.30
ps = {
"all p_i = 0.30": p_equal,
"half 0.10, half 0.50": p_mixture,
"half 0.00, half 0.60": p_extreme,
}
for name, pvec in ps.items():
m = poisson_binom_moments(pvec)
print(f"{name:>22} | mean={m['mean']:.2f}, var={m['var']:.2f}")
show_pmf_comparison(ps, title="Same mean, different heterogeneity → different variance/shape")
all p_i = 0.30 | mean=9.00, var=6.30
half 0.10, half 0.50 | mean=9.00, var=5.10
half 0.00, half 0.60 | mean=9.00, var=3.60
6) Derivations#
A) Expectation#
Because \(X = \sum_i X_i\),
B) Variance#
For independent \(X_i\),
C) Likelihood#
Suppose you observe i.i.d. samples \(x_1,\dots,x_m\) of the sum \(X\) (each sample is a sum of \(n\) Bernoulli trials with the same probability vector \(\mathbf p\)).
The likelihood is
Equivalently, the log-likelihood is
Two important practical notes:
The likelihood is invariant to permuting the entries of \(\mathbf p\) (only the multiset of probabilities matters).
Estimating a full length-\(n\) vector \(\mathbf p\) from only the sums is often ill-posed unless \(n\) is small or you add structure/priors.
def poisson_binom_log_likelihood(p: np.ndarray, data: np.ndarray) -> float:
p = validate_p_vector(p)
data = np.asarray(data)
n = p.size
if np.any((data < 0) | (data > n) | (data != data.astype(int))):
return -np.inf
_, pmf = poisson_binom_pmf_array(p)
logpmf = np.log(pmf)
return float(np.sum(logpmf[data.astype(int)]))
# A single-observation likelihood surface (n=2)
# This is small enough that we can visualize log L(p1,p2) on a grid.
x_obs = 1
p1_grid = np.linspace(0.01, 0.99, 80)
p2_grid = np.linspace(0.01, 0.99, 80)
LL = np.empty((p1_grid.size, p2_grid.size))
for i, p1 in enumerate(p1_grid):
for j, p2 in enumerate(p2_grid):
LL[i, j] = poisson_binom_log_likelihood(np.array([p1, p2]), np.array([x_obs]))
fig = go.Figure(
data=go.Heatmap(
x=p2_grid,
y=p1_grid,
z=LL,
colorscale="Viridis",
colorbar_title="log L",
)
)
fig.update_layout(
title=f"Log-likelihood surface for a single observation x={x_obs} (n=2)",
xaxis_title="p2",
yaxis_title="p1",
)
fig.show()
7) Sampling & Simulation (NumPy-only)#
The most direct sampler matches the definition:
For each trial \(i\), draw \(X_i \sim \mathrm{Bernoulli}(p_i)\).
Return \(X = \sum_i X_i\).
This is \(\mathcal O(n)\) work per sample and is typically fast in NumPy due to vectorization.
Below is a NumPy-only implementation using uniforms: \(X_i = \mathbb 1\{U_i < p_i\}\) with \(U_i \sim \mathrm{Uniform}(0,1)\).
def sample_poisson_binom_numpy(p: np.ndarray, size=1, *, rng: np.random.Generator) -> np.ndarray:
p = validate_p_vector(p)
n = p.size
if isinstance(size, (int, np.integer)):
size_tuple = (int(size),)
else:
size_tuple = tuple(size)
u = rng.random((*size_tuple, n))
return (u < p).sum(axis=-1)
p = np.array([0.10, 0.25, 0.60, 0.80])
s = sample_poisson_binom_numpy(p, size=10, rng=rng)
print("samples:", s)
samples: [3 2 2 2 1 1 3 2 1 2]
8) Visualization#
We’ll visualize:
the exact PMF (via dynamic programming),
the exact CDF,
and Monte Carlo samples from the NumPy-only sampler.
p = np.array([0.05, 0.10, 0.20, 0.35, 0.55, 0.70, 0.85])
k, pmf = poisson_binom_pmf_array(p)
cdf = np.cumsum(pmf)
mc = sample_poisson_binom_numpy(p, size=100_000, rng=rng)
emp = np.bincount(mc, minlength=p.size + 1) / mc.size
fig = make_subplots(rows=1, cols=2, subplot_titles=("PMF", "CDF"))
fig.add_trace(go.Bar(x=k, y=pmf, name="Exact PMF", opacity=0.65), row=1, col=1)
fig.add_trace(go.Scatter(x=k, y=emp, mode="markers", name="MC empirical PMF"), row=1, col=1)
fig.add_trace(go.Scatter(x=k, y=cdf, mode="lines+markers", name="Exact CDF"), row=1, col=2)
fig.update_layout(
title=f"Poisson binomial with n={p.size} (heterogeneous p_i)",
bargap=0.2,
)
fig.update_xaxes(title_text="k", row=1, col=1)
fig.update_yaxes(title_text="P(X=k)", row=1, col=1)
fig.update_xaxes(title_text="k", row=1, col=2)
fig.update_yaxes(title_text="P(X≤k)", row=1, col=2, range=[0, 1.02])
fig.show()
9) SciPy Integration#
SciPy (v1.15+) includes scipy.stats.poisson_binom.
Key methods:
pmf,logpmfcdf,sf(survival function; often numerically better for upper tails)rvs(sampling)stats(..., moments='mvsk')(mean/var/skew/excess kurtosis)entropy
We’ll compare SciPy’s PMF to our NumPy recursion and show typical usage.
About fitting: SciPy’s generic scipy.stats.fit does not currently support poisson_binom because the shape parameter p is vector-valued.
About fitting: poisson_binom has a vector-valued parameter p, and generic MLE fitting utilities are not currently designed for that case. In practice, you usually:
treat \(\mathbf p\) as known/estimated from other data, and then use the distribution for inference on the sum, or
impose structure on \(\mathbf p\) (e.g., a low-dimensional model) and fit that model.
p = np.array([0.05, 0.10, 0.20, 0.35, 0.55, 0.70, 0.85])
rv = stats.poisson_binom(p)
k = np.arange(p.size + 1)
pmf_scipy = rv.pmf(k)
_, pmf_numpy = poisson_binom_pmf_array(p)
print("max |pmf_scipy - pmf_numpy| =", float(np.max(np.abs(pmf_scipy - pmf_numpy))))
cdf_scipy = rv.cdf(k)
cdf_numpy = poisson_binom_cdf(k, p)
print("max |cdf_scipy - cdf_numpy| =", float(np.max(np.abs(cdf_scipy - cdf_numpy))))
mean, var, skew, excess_kurt = rv.stats(moments="mvsk")
print("SciPy mean/var/skew/excess_kurt:", float(mean), float(var), float(skew), float(excess_kurt))
print("NumPy mean/var/skew/excess_kurt:", poisson_binom_moments(p))
print("SciPy entropy (nats):", float(rv.entropy()))
# Sampling
r = rv.rvs(size=10_000, random_state=rng)
print("sample mean (SciPy rvs):", r.mean())
max |pmf_scipy - pmf_numpy| = 5.551115123125783e-17
max |cdf_scipy - cdf_numpy| = 1.1102230246251565e-16
SciPy mean/var/skew/excess_kurt: 2.8 1.1099999999999999 0.06926288077113398 -0.1184562941320082
NumPy mean/var/skew/excess_kurt: {'mean': 2.8, 'var': 1.1099999999999999, 'skew': 0.06926288077112831, 'excess_kurt': -0.1184562941319698}
SciPy entropy (nats): 1.4700953431353814
sample mean (SciPy rvs): 2.808
# Attempting to use scipy.stats.fit (as of SciPy 1.15)
#
# poisson_binom has a vector-valued shape parameter p, and SciPy's generic fitter
# isn't currently set up to optimize vector-valued shape parameters.
import warnings
p_true = np.array([0.2, 0.5, 0.7])
data = stats.poisson_binom(p_true).rvs(size=500, random_state=rng)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
try:
fit_res = stats.fit(stats.poisson_binom, data)
print(fit_res)
except Exception as e:
print("scipy.stats.fit is not currently supported for poisson_binom (vector-valued p).")
print("Error:", type(e).__name__, "-", e)
scipy.stats.fit is not currently supported for poisson_binom (vector-valued p).
Error: RuntimeError - The map-like callable must be of the form f(func, iterable), returning a sequence of numbers the same length as 'iterable'
# Optional: demonstrate a simple MLE when n is small (educational)
#
# Warning: estimating a full vector p from only sum data is often ill-posed.
# Here we fit a small n and show that recovery is up to permutation.
p_true = np.array([0.15, 0.35, 0.60, 0.85])
rv_true = stats.poisson_binom(p_true)
data = rv_true.rvs(size=2_000, random_state=rng)
def neg_log_likelihood_unconstrained(q: np.ndarray) -> float:
# unconstrained q in R^n -> p in (0,1) via logistic transform
p = expit(q)
ll = poisson_binom_log_likelihood(p, data)
return -ll if np.isfinite(ll) else 1e30
# initialize near a binomial approximation: all p_i equal to sample_mean / n
p0 = np.full_like(p_true, data.mean() / p_true.size)
q0 = np.log(p0) - np.log1p(-p0)
res = minimize(neg_log_likelihood_unconstrained, q0, method="BFGS")
p_hat = expit(res.x)
print("converged:", res.success)
print("true p (sorted):", np.sort(p_true))
print("mle p (sorted):", np.sort(p_hat))
print("note: order is not identifiable; only the multiset matters")
converged: True
true p (sorted): [0.15 0.35 0.6 0.85]
mle p (sorted): [0.4896 0.4896 0.4896 0.4896]
note: order is not identifiable; only the multiset matters
10) Statistical Use Cases#
A) Hypothesis testing (tail probabilities)#
If you know/assume the probabilities \(\mathbf p\), you can test whether an observed count \(x\) is unusually large/small.
Example: under the null model \(X \sim \mathrm{PB}(\mathbf p)\),
upper-tail p-value: \(\mathbb P(X \ge x) = \mathrm{sf}(x-1)\)
lower-tail p-value: \(\mathbb P(X \le x) = \mathrm{cdf}(x)\)
B) Bayesian modeling (uncertain probabilities)#
If each \(p_i\) is uncertain (e.g., you have a posterior over \(p_i\)), the posterior predictive distribution of the sum is a mixture of Poisson binomials.
A simple approach is Monte Carlo:
sample \(\mathbf p^{(s)}\) from the posterior,
compute the Poisson binomial PMF for each draw,
average the PMFs.
C) Generative modeling#
Poisson binomial is a natural building block for generative models that produce counts from heterogeneous binary events (e.g., conversions, failures, wins).
# A) Hypothesis test example: is an observed count unusually high?
p = np.array([0.02, 0.05, 0.08, 0.10, 0.15, 0.20, 0.25, 0.30])
rv = stats.poisson_binom(p)
x_obs = 5
p_value_upper = rv.sf(x_obs - 1)
print("Observed x =", x_obs)
print("Upper-tail p-value P(X >= x) =", float(p_value_upper))
Observed x = 5
Upper-tail p-value P(X >= x) = 0.00132964999999996
# B) Bayesian predictive example (Monte Carlo mixture)
#
# Suppose each p_i has a Beta prior and we observed (s_i successes out of m_i trials)
# in some historical data. We want the predictive distribution for the *next* round
# of n heterogeneous Bernoulli events.
from scipy.stats import beta
rng2 = np.random.default_rng(SEED)
n = 12
a0, b0 = 1.0, 1.0 # uniform Beta prior
m_i = rng2.integers(20, 80, size=n)
# synthetic historical success counts
p_latent = rng2.uniform(0.05, 0.8, size=n)
s_i = rng2.binomial(m_i, p_latent)
# posterior is Beta(a0+s_i, b0+m_i-s_i)
a_post = a0 + s_i
b_post = b0 + m_i - s_i
S = 5_000 # posterior draws
pmf_accum = np.zeros(n + 1)
for _ in range(S):
p_draw = beta.rvs(a_post, b_post, random_state=rng2)
_, pmf_draw = poisson_binom_pmf_array(p_draw)
pmf_accum += pmf_draw
pmf_pred = pmf_accum / S
k = np.arange(n + 1)
fig = go.Figure(go.Bar(x=k, y=pmf_pred))
fig.update_layout(
title="Posterior predictive distribution of the sum (mixture of Poisson binomials)",
xaxis_title="k (successes)",
yaxis_title="Predictive probability",
)
fig.show()
11) Pitfalls#
Invalid parameters: each \(p_i\) must lie in \([0,1]\). Values very close to 0/1 can make the distribution nearly degenerate.
Dependence: Poisson binomial assumes independent trials. Positive correlation typically inflates variance relative to the model.
Computing the PMF:
Naively summing over subsets is \(\mathcal O(2^n)\) and infeasible.
The dynamic-programming recursion is \(\mathcal O(n^2)\) and is fine for moderate \(n\).
For large \(n\), FFT-based methods can be much faster (SciPy uses fast algorithms internally).
Tail probabilities: prefer
sfover1-cdffor upper tails to reduce catastrophic cancellation.Fitting: estimating a full vector \(\mathbf p\) from only sum observations is typically ill-posed; consider adding structure (GLM, hierarchical priors) or using additional data.
Serialization: SciPy’s
poisson_binominstances currently do not support pickling/unpickling.
12) Summary#
poisson_binommodels the number of successes in independent, non-identical Bernoulli trials.The PMF is the coefficient of a polynomial / PGF; dynamic programming computes it in \(\mathcal O(n^2)\) time.
Mean and variance are simple sums: \(\sum p_i\) and \(\sum p_i(1-p_i)\); higher cumulants add as well.
NumPy-only simulation is straightforward: sample each Bernoulli and sum.
SciPy’s
scipy.stats.poisson_binomprovides PMF/CDF/SF/RVS/stats/entropy and is the go-to tool for inference on the sum.